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ABSTRACT 


This  report  documents  an  investigation  of 
three  recent  programs  for  the  solution  of  AX  = B 
where  A is  symmetric  and  sparse:  the  Yale  Sparse 
Matrix  Package;  the  Munksgaard  subroutines;  and 
the  Mesztenyi-Rheinboldt  subroutines.  The  first 
two  programs  compute  in-core  solutions;  the  third 
uses  random  access  to  obtain  its  solution.  The 
performance  of  these  three  programs  is  compared 
with  that  of  a previously  developed  out-of-core 
equation  solver,  CSKYDG2 . The  two  in-core 
equation  solvers  (especially  the  first)  are 
faster  than  CSKYDG2 ; the  third  is  slower.  All 
three  provide  the  same  degree  of  accuracy  as 
CSKYDG2 , but  they  require  large  amounts  of  core 
storage.  It  would  appear  that  the  two  in-core 
equation  solvers  are  not  suited  for  use  on  the 
CDC  6000  series  of  computers  in  their  present 
form  due  to  the  limited  amount  of  core  storage 
available.  Although  the  third  equation  solver 
does  not  require  as  much  core  storage,  it  does 
not  perform  as  well  as  existing  out-of-core 
equation  solvers  which  use  the  same  or  lesser 
amounts  of  core  storage. 


INTRODUCTION 


One  of  the  long  range  projects  of  the  Computation, 
Mathematics,  and  Logistics  Department  has  been  the  development 
of  mathematical  subroutines  suitable  for  use  in  the  computer- 
aided  structural  analysis  of  ships.  Many  unrelated  efforts 
in  both  government  and  industry  have  resulted  in  computer 
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programs  that  treat  particular  classes  of  structural  problems. 
These  programs  often  involve  the  solution  of  similar  mathe- 
matical problems  but,  since  the  solutions  are  reached 
independently,  the  efficiency  and  accuracy  of  the  various 
algorithms  used  may  vary  greatly.  The  need  to  coordinate 
these  diverse  efforts,  to  develop  improved  methods  of  more 
general  applicability,  and  to  produce  more  comprehensive 
programs  for  solving  Navy  structural  problems  became  obvious. 

A project  was  therefore  established  to  coordinate  research 
efforts  involving  mathematical  and  computational  methods  in 
the  area  of  structural  mechanics  and  to  integrate  the  work  of 
mathematicians,  computer  specialists,  and  structural  engineers 
in  this  field. 

The  present  considerable  interest  in  the  finite  element 
approach  to  structural  analysis  is  evidenced  by  the  wide- 
spread use  of  NASTRAN  (NAsa  STRuctural  ANalysis  program)  and 

other  such  programs.  According  to  the  NASTRAN  Theoretical 
1* 

Manual,  "From  a theoretical  viewpoint,  the  formulation  of  a 
static  structural  problem  for  solution  by  the  displacement 
method  is  completely  described  by  the  matrix  eguation  KU=P." 
Thus,  there  is  a need  for  accurate  efficient  computer  sub- 
routines capable  of  solving  these  large  sparse  positive 
definite  systems  of  simultaneous  linear  equations.  However, 
the  order  of  K is  often  so  large  that,  even  when  advantage 
is  taken  of  K's  symmetry  and  banded  structure,  it  is  not 
feasible  and  sometimes  not  even  possible,  to  store  K in  the 
core  memory  of  a computer . 

This  report  compares  three  recently  developed  programs 
for  solving  KU=P,  where  K is  a matrix  of  very  large  order, 
with  one  of  the  author's  previously  developed  out-of-core 
equation  solvers. 


* 


A complete  listing  of  references  is  given  on  page  35. 


THE  PROGRAMS 


The  following  programs  are  considered  in  this  report: 

2 3 

• The  Yale  Sparse  Matrix  Package  ' is  a collection 

of  subroutines  which  can  be  used  to  solve  both  symmetric  and 

non-symmetric  systems.  It  uses  algorithms  developed  by 

2 3 

Eisenstat,  et  al.  ' This  package  is  to  some  extent  proprie- 
tary since  limitations  are  placed  on  its  use  and  distribution, 

4 

• The  Munksgaard  collection  of  subroutines 
implements  Duff's  algorithm  for  solving  sparse  symmeyric 
systems.  It  can  obtain  a stable  decomposition  of  K in  both 

5 

the  definite  and  indefinite  cases.  This  program  was 
developed  at  the  Technical  University  of  Denmark  and  the 
Atomic  Energy  Research  Establishment  of  the  United  Kingdom. 

6 7 

• The  Mesztenyi-Rheinboldt  package  of  subroutines  ' 
can  be  used  to  solve  both  symmetric  and  non-symmetric  systems 
by  means  of  triangular  decomposition.  These  subroutines  are 
intended  for  use  with  systems  whose  coefficient  matrix  K will 
fit  into  core  but  may  not  do  so  after  decomposition.  Some 

of  these  subroutines  were  modified  by  the  present  author  to 
more  efficiently  utilize  the  random  access  storage  capabili- 
ties of  the  CDC  6000  series  of  computers. 

• CSKYDG2  is  an  out-of-core  Cholesky  algorithm 

g 

equation  solver  developed  by  the  present  author.  It  makes 
use  of  auxiliary  storage  via  the  random  access  capabilities 
of  the  CDC  6000  series  of  computers. 


THE  TEST  EXAMPLES 


EXAMPLE  1 


The  matrix  family  of  Table  1,  A^2  > generated  as 
follows:  Let  N be  an  integer  - 3.  Let  be  the  tridiagonal 
of  order  N with  4's  on  the  diagonal  and  a line  of  -I's  above 
and  below  the  diagonal.  Let  be  the  identity  matrix  of 
order  N.  An  (N+1) -banded  matrix  of  order  , A^  is  constructed 
by 

1.  stringing  N submatrices  along  the  diagonal, 

2.  inserting  lines  of  N-1  submatrices  above 

and  below  the  diagonal,  and 

3.  setting  the  remaining  elements  of  A^2  equal  to  0. 

Application  of  Gerschgor in ' s theorem  shows  AA2  to  be  positive 

^ 1 

definite.  The  right-hand  side  of  the  system  A^a^  = ^^^2  is 
chosen  such  that  all  components  of  the  exact  solution  vector 
except  the  first,  which  is  1,  have  the  value  0. 


EXAMPLE  2 

2 

The  matrix  family  of  Table  2,  A^2  » generated  as 

follows:  Let  N be  an  integer  - 3.  Assume  the  real  symmetric 

matrix  of  order  and  bandwidth  N with  on  the  diagonal 

and  -1  elements  filling  out  the  rest  of  the  band.  Then  change 

the  value  of  each  zero  element  in  the  last  N rows  and  columns 

2 

to  a -1.  As  before,  Gerschgorin ' s theorem  shows  A^2  to  be 

positive  definite.  Note  that  from  the  viewpoint  of  bandwidth, 
2 

Aj^2  is  ^ full  matrix.  The  right-hand  side  of  the  system 
2 

Aj^2^  " Djyj2  is  chosen  such  that  all  the  components  of  the  exact 
solution  vector,  except  the  f'  , which  is  1,  are  zero. 
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TABLE  NOTATION 


1 2 

The  solution  times  for  A.X  = and  A,X  = D are 

N N N N 

tabulated  in  Tables  1 and  2,  respectively.  The  CSKYDG2  times 

g 

were  originally  run  on  the  CDC  6600  at  DTNSRDC.  The  other 
programs  were  run  on  the  CDC  6400  at  DTNSRDC  and  the  resulting 
times  were  divided  by  3 to  give  the  equivalent  CDC  6600  times. 
The  column  headings  are  defined  as  follows: 

N - the  order  of  the  system 
M - system  bandwidth 

T - the  time  required  for  the  ORDV  subroutine  from 

the  Yale  package  to  reorder  the  system 
Tsolve  ~ time  required  for  either  the  SDRV  subroutine 

from  the  Yale  package  or  CSKYDG2  to  factor  and 
solve  the  system 

Ttotal  ~ total  time  required  by  a program  to  solve 

a system 

- the  time  required  for  either  the  INDANL 

DtiCUMF 

subroutine  of  the  Munksgaard  program  or  the 
SDECOl  subroutine  of  the  Mesztenyi-Rheinboldt 
program  to  factor  the  coefficient  matrix 
Tracksup  ~ time  required  for  either  the  INDOPR 

subroutine  of  the  Munksgaard  program  or  the 
SSLV  subroutine  of  the  Mesztenyi-Rheinboldt 
program  to  solve  the  factored  system 
Topmnn  " the  time  required  by  the  SETUP  preprocessor 

bfi T Ulr 

subroutine  for  CSKYDG2 


SOLVE 


TOTAL 


DECOMP 


BACKSUP 


SETUP 


(The  above  times  are  given  in  terms  of  CDC  6600  CPU  seconds.) 


TABLE  1 - EXECUTION  TIMES  FOR  THE  FIRST  SET  OF  SYSTEMS 
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OBSERVATIONS  AND  CONCLUSIONS 


An  examination  of  the  times  in  Tables  1 and  2 indicates 
that  the  Yale  program  is  the  fastest  by  far  of  all  four 
programs,  ranging  from  several  times  up  to  15  times  faster 
than  CSKYDG2  (total  time) . The  Munksgaard  program  is  at  best 
some  4 to  5 times  faster  than  CSKYDG2  (total  time)  for 
example  1 but  just  noticeably  faster  than  CSKYDG2  (total  time) 
for  example  2.  Why  the  Yale  program  is  so  much  faster  than 
the  Munksgaard  program  is  not  immediately  clear. 

The  Mesztenyi-Rheinboldt  program  appears  to  be  somewhat 
slower  than  CSKYDG2  because  CSKYDG2  moves  several  rows  at 
one  time  in  out  of  auxiliary  core  storage  whereas  the  present 
version  of  the  Mesztenyi-Rheinboldt  program  moves  only  one 
row  at  a time.  (As  a matter  of  fact  the  UNIVAC  version  of  the 
Mesztenyi-Rheinboldt  program  moves  individual  elements.  This 
was  changed  to  row  movement  by  the  present  author  when  con- 
verting the  program  to  the  CDC  6000  series  of  computers.) 

The  Yale  and  Munksgaard  programs  have  a tremendous 
time  advantage  over  the  Mesztenyi-Rheinboldt  program  and 
CSKYDG2  because  the  first  two  programs  do  everything  in  core 
while  the  latter  two  programs  use  random  access  to  move  rows 
in  and  out  of  auxiliary  storage.  However,  a field  length  of 
300000  CM  was  required  to  obtain  the  times  given  for  the  Yale 
and  Munksgaard  programs  in  Tables  1 and  2 even  though  the 
simplest  of  driving  programs  was  used.  (The  Mesztenyi- 
Rheinboldt  times  required  a field  length  of  225000  CM.) 

Clearly  the  use  of  some  device  such  as  "overlay"  would  be 
required  to  make  use  of  either  the  Yale  or  the  Munksgaard 
program  in  some  applications.  The  Yale  and  Munksgaard  pro- 
grams are  really  suited  for  computers  such  as  the  Texas 
Instruments'  Advanced  Scientific  Computer  which  have  abundant 
core  storage,  although  the  programs  would  have  to  be  rewritten 
for  the  most  part  to  obtain  the  full  benefit  of  the 
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computer's  optimizing  capability.  For  further  testing  on  the 
CDC  6000  series  these  programs  should  definitely  be  modified 
by  the  introduction  of  integer  packing  and  unpacking  sub- 
routines to  save  core  storage  by  cutting  dovm  the  size  of  the 
integer  arrays. 
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PROGRAM  LISTINGS 


Listings  are  provided  for  the  Munksgaard  and  Mesztenyi- 

Rheinboldt  programs  since  it  was  necessary  to  obtain  the 

times  in  Tables  1 and  2.  Listings  for  CSKYDG2  and  the  Yale 

8 2 

program  will  be  found  in  Gignac  and  Eisenstat  et  al. , 
respectively.  As  noted  previously,  the  Yale  Sparse  Matrix 
Package  is  considered  proprietary  in  some  sense. 


THE  MESZTENYI-RHEINBOLDT  PROGRAM 


SUMOUTINE  SI  NT  01 

• CONNON  NOIt)(FO(7)»RY(100(l0ltCT(10000)»O(10000)»AN(1600l  »NO(1600>  . 

COMMON  N0(0>tF0in  »RY I 23011) ^CV (23 00 0) I2i0 00) »*N(1600I >NO*i600) » 
1 1E(1600>«1H(1600) »IP(1600I «IN0EX( lEOl) *8 (1600 ) » X(160 0) 

INFEGER  RV,Cy 

C • INITIALIZE  SYMMETRIC  STRUCTURE  ANO  COEFFICIENT  ARRAYS  • 

C 

c 

N«RD(1) 

NO(S)«N 

F0(3)«0. 

00  10  IsltN 
AN(I)>0. 

RYCDsI 

CV(I)«I 

10  NO<I)3l 

RETURN 

ENB 


C 

C 

C 

c 


SUBROUTINE  SBLOOl (ItUt V) 

CONNON  NOIO)fFO(7)«RY(  10000)»CT(10000)  tA(10000)>AN(1600)  (NDdSOO)  * 
CONNON  NO(0)>FO(7)«RY(23000),CY(23  000> .A(23000) t AN(1600) (NO  11600 ) » 
1 IE(1600) t INI  1600) y IP (1600) tINOEX(1601) y8 (16 00) » X( 160 0) 

INTEGER  RTyCr 

* 8UILD  SYMMETRIC  STRUCTURE  ANO  COEFFICIENT  ARRAYS  • 


ANAX1(F0(3)  «ABS(  V)  ) 
AN(I)»S 

TO  20 


FO(3) 

AN(I) 

IF  (I.EQ.J)  GO 
N0(6)>M0(5)»1 
N0«ND(S) 

A(N0)>V 

AN( J)=AN( J)»S 

NO(I)«NO(I)»l 

NO(J)>NO(J)»1 

IF  (I.GT.J)  GO  TO  10 

RY(N0)3RV(I) 

CY(N0>«CYIJ) 

RY(I)3H0 

CT(J)>M0 

RETURN 


10  RY(N0)>RY(U) 
CY(N0)3Cy (I) 
RYI J)«M0 
CY(I)«N0 
RETURN 
20  A(I)«V 
RETURN 
ENB 
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noon  noo  non#  o o o o 


SUBROUTINE  SOECOi 

CONNON  NOIB)tFO(T)tRyi23000) tCV(230BO)tB(23000) «AN(1600I «NOC1600) » 
t IE(1600»*1H41600>»IPC1600) »1NDEXI 16011 tB(1600> «X< 1610 ) 

INTEGER  RV,Cy 
OINENSION  88(1600) 

CONNON  /ONE/  LIN1(1600> 


* OECONPOSE  SVHNETRIC  NATRIX  • 


N>N0( 1) 

NO(6>sO 

FO(S)*0. 

P0(6>>1. 

FO(2)«0. 

FO(6)>0. 

NNsO 

NO(6lsHO(5l 
NO(7) 3H0(S) 

N0(8l«0 
00  10  lsl»N 

FD(7)>F0(7)»At06( AN(I)  ) 

10  IP(I>-I 

FO(7|xO.$»FO(F) 

CALL  SVMI 

LOOP  ON  PIVOTING 

DO  260  I«1«N 
KxO 

IF  (I.EO.N)  GO  TO  30 

SELECT  PIVOT  BY  NININAL  DEGREE 

NOX*N»l 

AXsO. 

OO  1$  J«I»N 
IX«tP(J) 

IS  AXxAHAXK  AX«AeS(A  (IX)  I ) 

AX>AX*F0(2I 
OO  20  JxIfN 
IX«IP(J) 

IF  (AeS(A(IX)I.Lr.AX)  GO  TO  20 
IF  (NO(IX) .GE.NOX)  GO  TO  20 
NOKxNOdX) 

IY>J 

20  CONTINUE 

IF  (I.EO.IT)  GO  TO  30 
J«IP(IY) 

IP(iy)xIP(I) 

IP(I)xJ 

COLLECT  THE  ROM  AND  COLUMN  OF  THE  PIVOT 
ALSO  DELETE  THEN  FROM  THE  STORAGE 

30  IX>IP(1) 

S>A(1X) 

BBOVflMlE-COPt 
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o o o 


IF  c«es(si  .Lr.FOdii  go  to  sgo 

IF  lI.Ea.N)  GO  TO  Ul 

IYi«0 

IV>IX 

GO  IV«RV(1V) 

IF  <IV.EO.IX)  GO  TO  70 
IZ«1T 

50  IZ>CV(IZ) 

IF  (IZ.GT.NI  GO  TO  60 

KsK»l 

IE<K)>IY 

IHIX»*IZ 

ANIKMAdV) 

A(IT)«0. 

NOCIZ>>NO(IZI*i 
60  IF  fCYlIZi .NE.IVI  GO  TO  SO 
CV(IZ)«CY(IY) 

CY<IYI«NN 
NN>IY 
GO  TO  GO 
70  1Y*CTIIY> 

IF  (1Y.EQ.IXI  GO  TO  100 
IZ«IY 
lYt«IY 
SO  IZ«RT(IZ) 

IF  (IZ.GT.NI  GO  TO  90 

K«K»1 

IE(KI>IY 

IN(Kl3lZ 

AN(KI*A(1YI 

AdYlaO. 

N0(IZI=N0(IZ>>1 
90  IF  (RY(IZ) .NE.IY)  GO  TO  SO 
RYdZIsRYdYl 
GO  TO  70 

100  IF  (lYl.EQ.O)  GO  TO  110 
CYdYlMNN 
NNBCYdX) 


NOOIFICATION  OF  THE  RON  £I.EI(ENTS 


110  FO(GI«ANAXl(FO(GI  tABS(SI ) 
FO(5)»FO(5)»AI.OG(AaS(S» 

IF  (S.LT.O.I  F0(6l>-F0(6l 

CALL  SVH(-IX»St 

BB(1I«-IX 

BO(ZI*S 

LIN«G 

HO(S) 3MO(A|»l»K*K 
IF  (K.EQ.OI  GO  TO  250 
IF  (K.EQ.OI  GO  TO  2G9 
DO  115  JsltX 
AN(  JI>AN(  Jl/S 
call  SVM(IH(JI»AN(JII 
BB(LIN-1I>IH(  Jl 
BO«LIHI>AN (Jl 

FD(GI«ANAX1(F0(GI  »ABS(  AN(JIII 
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IIS  L1N*LIN»2 
Ki«K>l 

LOOP  FOR  THE  CROSS-POINT  ELENENTS 
DO  240  J>1»K 

Z««NIJ> 

A(IZI«A(IZ)-S*Z*Z 
F0(4)>ANAX1(F0(4) > ABSCACIZI )l 
IF  (J.EQ.K)  60  TO  240 
00  230 
JZ>IH(JJ> 

I1*NIN0(IZ»JZI 
I2*NAX0CIZ«JZ) 

Lt«RTCIll 
L2*Cy (I2> 

120  IF  KLl.EQ.Ill  .0R.(L2.EQ.12I  I GO  TO  140 
IF  (L1.EQ.L2)  60  TO  220 
IF  IL1.CT.L2I  60  TO  130 
L2«Cr(L2) 

60  TO  120 
130  L1«RV  (Lll 
60  TO  120 
C INSERTION  OF  A NON- ZERO  ELENENT 
140  NOf  I1)>‘ND(I1)  »1 
NOtI2l-NO(I2»  »1 
MO(7)«NO(TI«l 
IF  (NN.NE.OI  60  TO  170 
C USE  NEW  STORAGE 

NDt6)3N0(6)«l 

IF  (N0I6I  .GT.NOtZn  60  TO  310 
160  L1>N0(6) 

AdlMO 

RY(L1)3RV(I1) 

CVILllxCT(I2» 

RV(I1)3L1 
CYdZKLl 
60  TO  220 

C USE  AVAILABLE  STORAGE 
170  L1«HN 

HNsCr (HNI 

L3«I1 

L2<12 

180  IF  (RV(L3I  .LT.LII  60  TO  190 
L3>RV (L3t 
60  TO  ISO  ' 

190  RT(L1>3RTIL3) 

RV(L3)3L1 

200  IF  ICVtLZ) .LT.Ll)  60  TO  210 
L2sCT(L2l 
60  TO  200 
210  CV(L1»*CV(L2I 
CV(L2)>Ll 

C WRITE  OUT  CROSS-POINT  ELENENT 

220  AIL1>«AIL1»-ANIJ)«AN(JJ»*S  - — 

bki. 
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230  FOi<i)«ANAXtfFOUI  •ABSU(Lll)) 
2A0  CONTINUE 

C ENO  Of  NOOIFICATION  LOON 

• 250  CAUL  STM(-XX,SI 
2A9  BBILIN'IO-IX 

aB(LIN)>S 

LINI(1)«LIH 

250  CALL  NiarNSCF«B6*LIN»I) 

C 

C ENO  OF  PIVOTING  LOOP 
C 

• CALL  SOME 
RETURN 

C 

C SINGULAR  NATRIX 

C 

300  NO(l»)«l 
RETURN 
3i0  HOtk}s3 
RETURN 
C 

ENO 
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SUBROUTINE  SSLVTNOtXfV) 

DIMENSION  NO<t> «X(1) «V  (1> 

DIMENSION  88(16001 
CONNON  /ONE/  LlHKieOO) 

C * BACKSUaSTITUTION  FOR  SYNNETRXC  OECONFOSCO  MATRIX  * 

c • ••«•••••»•*»«»»«•»«»»••»•«•«««»•«•••«««••«••»«»«»» 

c 

N«iO(ll 
00  10  I«t.N 
10  V(I)»X(I) 

• CALL  SVRI 
I«0 

C FORNARO  SUBSTITUTION 

• 20  CALL  SMkFtL2$Zt 
20  I«I»l 

CALL  REAONSCTtBBtLINKlI  fll 

• L2*-L2 
L2*-BB(1I 
S«V(L2) 

LIN«<» 

• 30  CALL  StfRF(L2«Zt 
30  L2>BB(LIH>lt 

Z«BB(LINI 

IF  (L2.LT.0>  GO  TO  60 
r (L2»=»(L2»-Z*S 
LZNsLIN«2 
GO  TO  30 

40  IF  (I.LT.N)  GO  TO  20 
8ACKMAA0  BACKSUBSTITUTION 
SO  CALL  SVRB(L2»Z> 

• J»-L2 

SO  CALL  REAONS(/*Sa»LINI(  I»»I) 

LIN«LINI(II 

J»BB(L1M>1) 

Z*BB(L1N> 

\f(J»»VIJ»/Z 

I«I-l 

» 60  CALL  SVRa(L2»ZI 

60  LlN«LlH-2 
L2>BB(LIN>1) 

2«B8(LIN) 

IF  IL2.lt. 01  GO  TO  TO 
Y(JI*y<JI-Z*T<L2» 

GO  TO  60 

TO  IF  (I.GT.O)  GO  TO  SO 
RETURN 
C 

EN8 
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THE  MUNKSGAARD  PROGRAM 


SUBROUTINE  INDANLU(lN01»IN02tNZl>I  A1»I  A2*N»1N»IP»0»C»U«NI 
INPLICir  REAL'S  (A-H,S-Z) , INTEGER  (I-RI 
INPLICir  REAL  (A-t*,S-Z)  .INTEGER  (I*RI 
DOUBLE  PRECISION  A(IA2).0CN) . H(Nt 
OINENSION  A(IA2lf OCN)tN(NI 

INTEGER  IP(N,2) .KP(2> . KL (21 . JP(2) .FIL 
INTEGER'Z  IN(N»n  .INOl  (lAll  , IN02  (IA2I 
INTEGER  INCN.TI ,IN01  ilAl) .IN02<IA2I 
LOGICAL  FIRST, STA81.STAB2 
CONNON/IN  OEFI /F IL . NUEL .LROM.LCOL . NCP ,NUCL 
IPd.ll.IPd.Z)  POINT  TO  THE  START  OF  RON/C3LUHN  I. 

IN(I.1I,IM(I.21  HOLD  THE  NUNBER  OF  NONZEROS  IN  ROH/COLUNN  I OF  THE 
LONER  TRIANGULAR  PART  OF  A. 

DURING  THE  NAIN  BOOT  OF  THIS  SUBROUTINE  THE  TECTORS  IN(',  3).  IH(«.<») 
INI'.SI  ARE  USED  TO  HOLD  OOUBLT  LINKED  LISTS  OF  RONS  THAT  HAVE 
NOT  BEEN  PIVOTAL  AND  HAVE  EQUAL  NUMBER  OF  NONZEROS. 

INd.3)  HOLD  FIRST  RON/COL UHN  TO  HAVE  I NONZEROS  OR  ZERO  IF  THERE 
ARE  NONE. 

INd,<»l  HOLD  RON/COLUHN  NUMBER  OF  RON/CBLUHN  PRIOR  TO  RON  I IN  ITS 
LIST  OR  ZERO  IF  NONE. 

INd.SI  HOLD  RONFCOLUHN  NUMBER  OF  RON/COLUNN  AFTER  RON  I IN  ITS  LIST 
OR  ZERO  IF  NONE. 

IN<*.6>  AND  INI', FI  ARE  USED  TO  UNPAC  THE  PIVOT  RON(S)  INVOLVED  IN 
AN  ACTUAL  RON  OPERATION. 

DURING  THE  MAIN  BOOT  OF  THE  SUBROUTINE  INOl  AND  IN02  KEEP  A COLUMN 
FILE  AND  A RON  FILE  CONTAINING  RESPECTIVELT  THE  RON  NUMBERS  OF 
THE  NON-ZEROS  OF  EACH  COLUMN  ANO  THE  COLUMN  NUMBERS  OF  THE  NON- 
ZEROS OF  EACH  RON.  THE  NON-ZEROS  OF  A FOLLONS  THE  ORDERING  OF  THE 
RON  FILE  IN02.  THE  IP  ARRAYS  POINTS  TO  THE  START  POSITION  IN  IN02 
(ANO  Al  ANO  INDl  OF  EACH  RON  ANO  COLUMN. 

00  5 Isl.N 
00  A J=6,Z 

4 IN(I,J)s-l 

0(1) -0. 

NdlsO. 

DO  5 J=1«S 

5 IN(I.J)=0 
GsQ. 

NZsNZi 

NUEL^NZ 

L = l 


COUNT  NUNBER  OF  ELEMENTS 
00  20  lOUNNYsl.NZ 
IF  (L.GT.NUEL)  goto  25 
00  10  KsL.NUEL 
1=1ND1(K) 

JsIN02(O 

IF  (I.LT.l  .OR.  I.GT.N)  GOTO  610 
IF  (J.LT.l  .OR.  J.GT.NI  GOTO  610 
» Gl30ABS(A(KII 

GIz  A8S(A(K>) 

* G^ONAXKGl.G) 

G>AMAX1(GI .G) 

IF  (I.EQ.J)  GOTO  15 
IF  (H(I).LT.GI) N(l)sGI 


18 


ouu  uuu  oou  uu 


IF  (M(J).LT.GI)M(J)«GI 

10  IMU,2)>lMIJ,2l*i 
GOTO  25 

CHECK  FOR  DOUBLE  ENTRIES  ON  THE  DIAGONAL  AND  REMOVE  DIAGONAL  FORM 
THE  file. 

15  IRsl 
J»I 

IF  (0(11 .NE..0I  GOTO  650 
0(II=A(K) 

L»K 

A(L>*A(NUEL» 

INOKLXINOKNUELT 

IN02(L):IND2(NUEL) 

IN01(NUEL>*0 
IN02(N(CLI«0 
20  NUEL«NUEL*1 

NCP  IS  THE  NUNBFR  OF  COMPRESSES  PERMITTED  BEFORE  A MORNING  RESIR.TS. 
NCP  IS  THE  NUMBER  OF  COMPRESSES. 

25  NCP*0 

NCP«NAX0(N/10«2a) 

NZ^NUEL 

LCOLsNUEL 

NUCLsNUEL 

LROM=NUEL 

CHECK  FOR  NULL  RON  ANO  INITIALIZE  IP(I»A»  AND  1P(1»2I  TO  POINT  JUST 
8EV0NT  NHERH  THE  LAST  COMPONENT  OF  RON/COLUMM  I OF  A MILL  BE  STORED. 
KI-1 

KJsi 

DO  30  I-1»N 
KIsKI»IM(I«l) 

KUsKJ>IM(I,2) 

IP(I»l|sKI 

IP(I*2)-KJ 

IF(IM(I,1I  »IM(I»2) .LE. 0 .AND.  0(1). £0.0.1  GOTO  630 
30  CONTINUE 


REORDER  BY  RONS  USING  IN-PLACE  SORT  ALGORITHM. 

00  50  I>1*NZ 
IRlsINOKI) 

C IF  IRl  IS  NEGATIVE  THE  ELEMENT  IS  IN  PLBCE  ALREADY. 
IF  (IRl.LT.OI  GOTO  50 
IC1>IN02(I) 

AlsA(I) 

Kl3lP(IRl>t)-l 
00  65  IOUNNV>l»NZ 
IF  (I.EQ.Kt)  GOTO  66 
lR2>IN01(Kt) 

IC2sIN02(Kl) 

A2sA(Kl) 

A(K1)3A1 
IND1(K1)3-IR1 
IND2(K1I«IC1 
IP(ZR1»1)«K1 


BtSLAVAllABUCOPY 


19 


1R1=IR2 

IClsICZ 

A1=A2 

•»S 

46  IF  (lOUMMY.EQ. 11  GOTO  49 
AlIlsAl 
1N01<I)=-IR1 
IND2a)  = ICl 

49  IPllRltDsI 

50  CONTINUE 
C 

C CHECK  FOK  OOUBBEL  ENTRIES  WHILE  USING  THE  CONSTRUCTED  ROM  FILE  TO  SET 
C UP  THE  COLUiN  FILE, 

KLLsNZ 
DO  60  I=1,H 
IR=N*1-I 
KPP=IPaR»l) 

IF  IKPP.GT.KLL)  GOTO  60 
00  55  KsKPPtKLL 
JsIN02(K) 

IF  (IN(J,4).EQ.IR)  GOTO  650 
IN(J,4I=IR 
KR=IP(J,2»-1 
IP(J,2>sKR 
55  IN01IKRI=IR 
60  KLL=KPP-1 
C 

C SET  UP  LINKED  LISTS  OF  ROHS/COLUMNS  WITH  EQUAL  NUMBER  OF  NON*ZEROS. 
DO  62  1=1, N 
N2I=IH(I,ll»IN(l,2t»l 
IN=  IN(NZI,3) 

IN<NZI,3i=l 

INII,5I=IN 

IN(I,4)=0 

62  IF  (IN.NE.O)  1N(IN,4I  = I 
C 

C START  THE  ELZHZNATI ONLOOP 
C 

PP=0 

GG=G 

DO  500  IIP=1,N 
IF(PP.NE.2)  GOTO  70 
PP  = 0 
GOTO  500 

70  MV1=<N-1»**2 

MY2=(2*N-4I*(N-2I 
MYM=MY2 
FIRSTS. TRUE. 

STAB2= .FALSE. 

C 

C SEARCH  RONS  WITH  RJPl  NONZEROS. 

DO  145  RJPl=i,N 
C 

C HYl  AND  HT2  ARE  THE  LEAST  COSTS  SO  FAR  FOR  STABLE  1X1  AND  2X2  PIVOTS. 
C HCH  IS  THE  THE  LEAST  OBTAINESLE  COST. 

HCM*MIN0( (RJP1-1»**2, I2*RJPl-4»**2/2» 

IF  CNYH.GT.MCH)  GOTO  72 
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IF  (Hyi.GT.MVZ. OR. FIRST)  GOTO  FI 
GOTO  210 

FI  IF  (.NOT.  STAB2>  GOTO  F2 
GOTO  220 

F2  JJP-IH  (RJP1,3I 

LOOP  ON  ROUS  OF  RJPl  NONZEROS. 

00  140  IDUNHT  sl,N 
IF(JJP.LE.O)  GOTO  14S 
111  = 0 

MVB=(2*N-4l**2/2 
IF  (.NOT.  FIRST)  GOTO  F5 
AU=0A8S(0( JJP) )/U 
AU=  AdS(0(JJP))/U 
STAB1=(H( JJPI .LE. AU) 

IF  (STAOl)  GOTO  105 
7S  P2=0 

COMPRESS  KOUFILE  IF  THERE  IS  NOT  ROOM  FOR  NEM  ELEMENTS. 

IF  (NUEL»RJP1.L T. IA2)  GOTO  F5 
IF  (LKOU»RJPl.GT. IA2.0R.NCP.GE.HCP)  GOTO  6F0 
CALL  C0MPRE(A,IND2>I A2  .N >IM ,I P » .TRUE . ) 

PERFORM  THE  1X1  STABILITY  TEST  OF  ROM/C«LUMN  JJP. 

F5  OO  86  L=lt2 
KPP=1P(JJP,L) 

KLL=KPP*IU(JJP,L) -1 
IF  (KPP.GT.KLLI  GOTO  S6 
OO  85  K=<PP,KLL 
IF  (L.EQ.Z)  GOTO  77 
PP  = IN02('() 

GOTO  78 
77  PP=IN01(K) 

F8  RJP2=IM(PP,1» ♦IU(PP,2) *1 

BUILD  UP  THE  FULL  RON  JJP  AT  THE  END  OF  THE  ROM  FILE.  OFF  DIAGONAL 
ELEMENTS  NHICH  HAS  BEEN  INVES TEG ATEN  AS  2X2  PIVOT  CANDIDATES  ARE 
SET  TO  THE  NEGATION  OF  THEIR  COLUMN  NUMBER. 

IN02(IA2-I11)  =-PP 
IF  (RJP2.lt. RJPl)  GOTO  85 
IF  (RJP2.GT.RJP1)  GOTO  83 
JLK=lH(JJPf4) 

OO  82  KOUHMT  =1,N 
IF  (JLK.EQ.O)  GOTO  83 
IF  (JLK.EU.PP)  GOTO  85 

82  JLK=IH(JLK«4) 

83  lN02(IA2-lIl)sPP 
MCN=(RJPl»RJP2-4) **2/2 

IF  (NCN.GE.MT2.0R.NCN.GE.NTB)  GOTO  85 

NTB=MCN 

P2  = PP 

KKP=K 

85  II1=II1»1 

86  CONTINUE 
MCN=NV2 

OO  104  JOUNNVsI.N 
IF  (P2.LE.0)  GOTO  110 
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CHECK  THE  STAaiLITV  OF  THE  2X2  PIVOT  OETERNINEO  BT  RON  JJP  AND  P2. 
JP(1)=JJP 
JP(2)=P2 
II2«IA2 
AAl^OCJPCin 
AA2=0(JPI2n 

IF  (IN02(KKP)  .EQ.P2)  GOTO  8« 

KPP*IP(P2,1» 

KPL=KPP*lN(P2,l»'l 
00  *2  KKP=<PP,KPl 
IF  (IN02(KKPI  .EQ. JP( 1) ) GOTO  S8 
82  CONTINUE 
88  AAt2=A(KKP) 

0ETP=AAl*AA2-AA12**2 
IF  (OETP.NE.O.  I GOTO  91 
APU-1.0020 
GOTO  92 

91  APUsOHAXKOABSI  AAll«0AaS(AA2l ) 

91  APU=ANAX1(  ABSCAADt  A8S(AA2M 
APU=<APU»0ABS(AA12)l*U/0AaS(0£TP) 

APU=(APU«  ABS<A A12)) *0/  ABS(OETP) 

DETERNINE  THE  NUMERICAL  MAXIMAL  ELEMENT  OF  RONS  JPdIfANO  JP(2I 
NHICH  IS  NOT  CONTAINED  IN  THE  BLOCK  DIAGONAL. 

GI  = 0MAX1(H  UP  III)  .HUP  (2))) 

GI^AHAXIIH  UP  (in  .NUP  (2)1) 

IF  ( APU’SI.GT.l)  GOTO  92 

JP1=JJP 

JP2=P2 

KP12=KKP 

A11=AA1/0ETP 

A22=AA2/0ETP 

A12=AA12/0ETP 

MY2  = Mt  B 

STAB2=.TRUE. 

GOTO  110 

92  PP2=0 
KB=0 
JPL=JJP 

00  98  L=1.2 
KPP=IP(JPL.L) 

KLL=KPP*IH(JPL.L) -1 
IF  (KPP.GT.KLL)  GOTO  98 
DO  92  K=KPP,KLL 
PP=IN02(II2) 

IF  (PP.EQ.JP(2) ) GOTO  95 

IF  (PP.LE.O)  GOTO  96 

HC=(RJP1*IH(PP,  1)  UNIPP,  2) -3)  **2/2 

IF  (MC.GE.NCN)  GOTO  96 

HCN=MC 

KB=K 

PP2=PP 

GOTO  96 

95  IND2(II2)*-P2 

96  112=112-1 
92  CONTINUE 


98  CONTINUE 
HVB^NCN 
NCNSNY2 
P2=PP2 
lOH  KKP^KB 
C 

C IF  THE  FIRST  STABLE  1X1  PIVOT  IS  FOUND  STORE  RELEVANT  POINTERS. 

105  NT1=(RJP1-1)**2 

JPP=JJP 
F1RST=, FALSE. 

C 

C OETERHINE  NETHER  THE  SPARSITY  CONDITION  IS  SATISFIED  OR  NOT. 

110  HYN=NIN0(HY1«NT2) 

IF  (NYN.GT.HCH)  GOTO  140 
IF  (NY1.GT.NY2. OR. FIRST!  GOTO  120 
GOTO  210 

120  IF  (.NOT.  STABE)  GOTO  140 
GOTO  220 

140  JJP-IM(JJP,5I 
145  CONTINUE 
C 

C PIVOT  IS  FOUND 
C 

C RON  JPl  IS  USED  AS  1X1  PIVOT. 

210  JPl=JPP 
PPsl 

GOTO  230 
C 

C RONS  JPl  ANO  JP2  ARE  USED  AS  2X2  PIVOT. 

220  PP=2 

JP(2)=JP2 
230  JP(1)=JP1 

C 

C REHOVE  RONS/COLUNNS  INVOLVED  IN  ELININATION  FROH  ORDERING  VECTORS. 
00  266  L1=1.PP 
JPL=JP(L1» 

DO  250  L=lt2 
KPPsIP(JPL,L) 

KLLslN(JPL.L)  ♦KPP-1 
IF  (KPP.GT.KLLI  GOTO  250 
00  246  K=KPPfKLL 
JSIN02 (K) 

IF  (L.EQ.2!  J>1N01(K) 

IL=IN(J,4I 
IN::IM(  J,5I 
IN(J,5>>-1 

IF  (IN.LT.OI  GOTO  246 
IF  (IL.E3.0)  GOTO  240 
lN(lLt5)-IN 
GOTO  245 

240  NZ=INIJ,1I*IN(J*2)»1 
IH(NZt3)sIN 

245  IF  (IN.GT.OI  IN(IN«4>«IL 

246  CONTINUE 
250  CONTINUE 

C 

C RE  HOVE  JPILII  FRON  ORDERING  VECTDRS 
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IL=1M( 

INsINIJPLfS) 

lH{JPL.9)s-l 

IF  IIN.LT.OI  COTO  266 

IF  (IL.EQ.O)  GOTO  260 

IH(IL,5I=IN 

GOTO  265 

260  NZ-IM(JPL,1)*IM(JPL»2) *1 
lN(NZ,3lslN 

265  IF(IN.GT.O) 

266  CONTINUE 
C 

C STORE  PIVOT. 

00  26T  L=l,PP 

267  INCJP(L) .4)=*11P»1-L 
C 

C CREATE  A NEM  ENTRY  FOR  RON  JPl  (AN0JP2i  IF  NESCCSSARV  ANO  TRANSFORN 
C COLUNN  PART  TO  ROMFILC.  RENOVE  PIVOTAL  ROV(S»  FROM  THE  COLUMN  FILE. 
DO  325  LL>1|PP 
JPL=JP(LL> 

NZCsIMCJPL  t2> 

IF  (NZC.EQ.OI  GOTO  290 
C 

C COMPRESS  RONFILE  IF  NECESSART 

IF  (NUEL«NZC*IMUPLtl)  .LT.IA2)  GOTO  270 
IF  (LR0H»NZC»IH(JPL«1> .GT.IA2  .OR.  NCP.GE.NCP)  GOTO  670 
CALL  C0NPRE(A,IN02yIA2rN»lM,IP».TRUE.) 

270  KPP=IPUPL,2) 

KLL=KPP*MZC-1 
IPIJPLy2>  sNUEL»l 
00  200  K=KPP,KLL 
NU£L3NUEL»1 
LC0L«LC0L>1 
IsINOKK) 

KIP>IP(I,1I 
KILsKIP^INd.l)  •! 

DO  275  K<3KIP»AIL 
IF  (IN02(KIC>.EQ.JPLI  GOTO  276 

275  CONTINUE 

276  1N02(NUEl)«I 
A(NUELI«AUIO 
INOKKIsO 

C 

C MOVE  ELEMENT  FROM  RONFILE. 

KRL=IP<l«l)»IHCItll-l 
C MOVE  ELEMENT  FROM  RONFILE. 

KRL*IP(Itl>»IM(I|ll*l 
IF  (KK.EQ.KRLI  GOTO  279 
ILslN02fKRLI 
IN02(KK|3lL 
A(KKi«A(KRL» 

279  IND2<KRL)sO 

200  INfI,l)«lNIIfl>*l 
C 

C MOVE  RONFILE  OF  ROM  JP(LL)  IF  NZC  .GT.  0 

290  KPP=IP(JPL«1I 

IF  (NZC.GT.OI  IPUPL«1>>IP(JPL»2I 
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KLLsiMUPL  ,1>^KPP-1 
IF  (KPP.GT.KLLl  GOTO  320 
DO  319  K:<PP>KI.L 
J=1N02(K) 

KPCsIP(J.2} 

NZ=IH(J,2>>1 

IM(J,2>=NZ 

KiCsKPC»N2 

IF  (KLC.GE.KPCI  GOTO  31J» 

1N01(KPC>=0 

GOTO  317 

314  00  319  KC=KPCt  KLC 

IFUPL.EO.lNDlIKCn  GOTO  316 
319  CONTINUE 

316  IN0l(KCt=lN01(KLC> 

INDKKLCl  =0 

317  LCOLsLCOL-1 

IF  (NZC.EQ.OI  GOTO  319 

NUEL3NUEL»1 

A(NUEL)=A(K) 

IN02(NUEL)-IN02(K) 

lND2(K>=a 

319  CONTINUE 

UPDATE  ROMPOINTE^S  TO  NOW  JPILL) 

320  IH(JPL,li»IN(JPC»l)«N2C 
329  IN( JPL,2J s-PP 


NOVE  THE  PUOrJkL  OFF  DIAGONAL  ELENENT  TO  THE  FRONT  OF  RON  JPl  IF  PP>2 
NZC=IM<JPlyl> 

IF  <PP.ea.ll  GOTO  360 
KPP=IP<JP1,1» 

KLL=KPP*NZC-1 
00  326  K=KPP><LL 
J-IN02(K> 

IF  iJ.Ea.JPZ)  GOTO  327 

326  CONTINUE 

327  IF  (K.EQ.KPPt 
Al-A(KPP) 

IlsIN02(KPP) 

A(KPP)=A(K) 

IN02(KPPI 3JP2 
ACKMAl 
IN02<K>sIl 


GOTO  32B 


IF  PP32  THEN  CONSTRUCT  A PSEUOORON  CONTAINING  THE  UNION  OF  COCUNN 
NUH8ERS  OF  THE  TNO  PIVOTAL  PONS  AT  THE  END  OF  THE  RONFILE* 


COMPRESS  RONFILE  IF  THERE  IS  NOT  ROOM  FOR  NEN  ELEMENTS. 

328  NZC3lW(JPl,l)*ZN<UP2«17'l 

IF  (NUEL«NZC.LT.IA2)  GOTO  329 
IF  <LR0N*NZC.GT.  IA2  .OR.  NCP.6E.NCP>  GOTO  670 
CALL  C0NPRE(A,INO2«IA2«NtINtIP..TRUE.) 

329  KPl3lP{JPl,l)»i 
KL1=IN(JP1,1>  »KPl-2 

IF  (KPI.GT.KLII  GOTO  331 


00  330  K^KPltALl 
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330  INfIN02(K> tO}*! 

331  NZC*IA2*1 
KPZ=IPCJP2,1I 
KL2«IM(JP2ttl  »KP2-1 

IF  IKP2.Gr.KL2l  GOTO  341 
00  340  K:KP2,KL2 
l3lN02(<> 

MZC^MZC-l 

IN02INZC1=I 

340  1M(I,6>>-1 

341  IF  (KPI.GT.KLII  GOTO  3S1 
DO  350  Ks<Pl,KLl 
IsIN02(KI 

IF  llMfIt6I.EQ.-ll  GOTO  350 

NZC<NZC-1 

IN02INZC)«I 

349  IMIIt6)--l 

350  CONTINUE 

351  NZCsI A2-NZC*! 

PERFORN  THE  ELIMINATION  LOOPING  ON  NONZEROS  IN  PIVOT  COLUNNTSl. 

360  IF  (NZC.EO.O)  GOTO  405 
DO  363  L=ltPP 
KPIDsIPlUPIDt  1) 

KL(LI  = KP(L)^1M(  JP(L)  «1)-1 
IF  (PP.EQ.Z.ANO.L .EQ.l I KP(L)«KP(LI*1 
IF  (KP(L)  .GT.KL  (LI  I GOTO  363 

UNPACK  PIVOT  RON(SI  IN  IN(»y5*LI. 

KPP=KP(LI 

KLLsKL(L) 

H(JP(L)>s0. 

DO  362  K=KPPtKLL 
J=IN02(K) 

N(JI-0. 

362  IM(J,5»L) sK-<PP 

363  CONTINUE 

DO  400  NC^ltNZC 
KC=IP(  JPlt  1)»NC-1 
IF  (PP.E0.2)  KC=IA2-NZC»NC 
IRslN02(KCI 

IF  <PP.E3.1I  AL3A(KCt/0( JPl) 

COMPRESS  RONFILE  UNLESS  IT  IS  CERTAIN  THAT  THERE  IS  ROOM  FO^  NEH  RON. 
IF  (NUEL»IN(IR, 1) ♦NZC.LE.IA2-N2C)  GOTO  365 
IF  (NCP.CE.NCP  .OR.  LRON«lN(IR> II i-NZC.GT . lAZ-NZCI  GOTO  620 
NZ0«NZC 

IF  (PP.EQ.ll  NZ0«0 

CALL  COMPRE(A,IN02.IA2-NZO.N,IN,IP..TRUE.) 

00  364  LslfPP 
KPILI^IPI JP(Llt It 

364  KL(LI«KP(LI»IN(UP(LI tll-1 
IF  (PP.E0.2I  KP(1I>KP( II ♦! 

365  KR<IP(IR«1I 
KRL«KR»IN(IRtlt-l 
IP(IR,1I«NUEL»1 

IF  (PP.EQ.ll  GOTO  372 
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AIRl-0. 

AIR2>0. 

IN123lM(IR,7) 

IF  (IMIl.GE.OI  AIRlsAf  KPtD^IMItl 
IF  CIMI2.GE.0I  AlR2«AIKPf2l«lMl2) 


TRANSFER  NOOIFIEO  ELEMENTS. 

77  IF  (KR.GT.KRLI  GOTO  420 
DO  390  KS=KR,<RL 
JslND2(KSt 
AJRlxO . 

INJl^INf J,6> 

IM(J»6>=-IMJl-2 

IF  (IMJ1.6E.0I  AJRlzA(KP(tl*IMJl» 

IF  (PP.E3.1)  GOTO  3S0 
AJR2=0. 

IMJ2=IHIJ,7) 

IMlJ,7»=-IMJ2-2 

IF  (IMJ2.GE.0I  AJR23A(RP(2I»IMJ2) 

A1  = A(KS)-AIR1*AJR1*A22UIR2*AJRI*A12-AIR2*AJR2*A11»AIRI*AJR2*A12 
GOTO  3S1 

ISO  A1  = A<KS>-AL»AJR1 
3S1  INO21KS)s0 
GI=0ABS(A1I 
GI-  ABS(Al) 

IF  (H(IR)  .LT.GII  HdRIsGI 
IF  (M(J). LT.GII  MtJIsGI 
G=0NAX1(G«GI) 

GsANAXl(G,GII 
NUEL=NUEL«1 
AtNUELIsAl 
390  1N02(NUELI=J 


SCAN  PIVOT  RON  FOR  FILLS. 

420  00  480  NCl3l,N2C 

IF  (PP.EQ.l)  GOTO  421 

KS=IA2*1-NC1 

GOTO  422 

421  KS»KP(1) -1*NC1 

422  J=IN02(KSI 
AJRlsQ. 
lMJlsIH(J,6) 

AJR2sQ. 

INJ2=IH(J,7) 

IF  <1MJ1.LE.-1  .AND.  1MJ2.LE.-1  .OR.  J.GT.IRI  GOTO  469 
IF  (INJl.GE.OI  AJR1«A(I(P(1I»INJ1> 

IF  (PP.EQ.p  GOTO  42S 
IF  (IMJ2.GE.0)  AJR2«A(KP(2)»IMJ2I 

41=-A1R1*AJR1*A22*AIR2*AJR1*A12-AIR2*AJR2*611*AIR1*AJR2*A12 


GOTO  42b 

425  A1«-AL*AJR1 

426  IF  (IR.NE.J)  GOTO  429 
0(IR)«0(IRI*A1 
AlsO(IR) 

• GsONAXltOABSIAllyGI 

GsANAXK  ASStAlltG) 
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Goto  G69 

G29  NUEL«NUEL»1 
LROM«LRON»1 
A(NUEL)>«1 
IN02(NUELI«J 
CIxOABSUll 
61*  A8S<«1> 

IF  (N(IR)  .LT.CI) 

IF  (M(J».LT.GII  MiJIsGI 

CREATE  FILL  IN  COLUHN  FILE. 

N2«IMU,2) 

KslPCJtE) 

KLt>K*NZ-l 

IF  POSSIBLE  PLACE  NEH  ELEMENT  AT  THE  END  OF  PRESENT  ENTRY. 
IF  TKLl.NE.NUCLI  GOTO  A30 
IF  ILCQL.GE.IA2I  GOTO  <«4|0 
NUCL*NUCL»1 
GOTO  G39 

L30  IF  (INOKKLl^lT  .NE.B)  GOTO  AAA 
l»35  IN01(KL1»1»>IR 

LCOL>LODL»t 
GOTO  AGS 

NEN  ENTRY  HAS  TO  BE  CREATED. 

|>A0  IF  INUCL»N2*l.Lr.ZAl)  GOTO  ASO 

COMPRESS  COLUMNFILE  IF  THERE  IS  NOT  ROOM  FOR  MEM  ENTRY. 

IF  (NCP.GE.NCP  .OR.  LCX»NZ«1  .GE.IAil  GOTO  GB2 
CALL  C0NPREIA»IN01»IA1,N,1M(1»2I .IPt 1»2I » .FALSE. ) 
K«IP(3,2> 

KLisK«NZ>l 

TRANSFER  OLD  ENTRY  INTO  MEM. 

^SA  lP(Jf2)sNUCL»l 

IF  (K.GT.KLII  GOTO  A61 
OO  AGO  KKsKfKLl 
NUCL>NUCL»1 
INOKNUCLKINOi  (KKt 
AGO  INOKKKTsO 

AOO  THE  NEH  ELEMENT. 

A61  NUCL-NUCLAl 
IN01<NWL»«IR 
LCOLsLCOL»t 
AGS  G«DNAX1IG«G1> 

AGS  G>AHAX1IG»GII 
INI J|2)«NZ»L 

AG9  if  fIM(J»ri.Lr.«lt  IM< J,7)«-lH(J»r}-2 
IF  (IN(J,6).LT.-1)  IN(J,G)>-1M(J,6)»2 
lMfIR»l)>NUEi.*I-lP(IR»l) 

AAA  CONTINUE 

CLEAN  IM<*,G)  ANO  INI«»7I  BY  SCANNING  TNE  PIVOT  RONS 
OO  AAA  NC«t»NZC 
IFIPP.EQ.ll  GOTO  AAl 
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K=IA2»1-NC 

GOTO  482 

481 

K=KP(1)*MC-1 

A(K)=A(<)/0(JP1) 

482 

JsIND2(K) 

00  484  L=lyPP 

484 

IM(J,5»L)=-1 

C INSERT  ROMVCOLUriNS  INVOLVEO  IN  EcIHlNATION  IN  LINKED  LISTS  OF  EQUAL 
C NUMBER  OF  NON  ZEROS 
485  Kt=IP(JPl,t) 

IF  (PP.EQ.Z)  <l=IA2-NZC*l 
K2=K1*NZC-1 

IF  (K1.GT.K2I  GOTO  A91 
00  490  K=K1,K2 
IK=IN02(K) 

NZ=IN(IR,1I»IH(IR.2I»1 

IN=IH(NZi3) 

IM(IR,5I=IN 

1NIIR,4)=0 

IH(NZ.3>=IR 

490  IF  (IN.NE.O)  IK(IN,4)  = IR 

491  IF  (PP.EQ.l)  GOTO  500 
0(JP1)=A22 

0( JP2)=A11 
4(IP<JP1,1H=-A12 
500  CONTINUE 

C ELIMINATION  LOOP  TERMINATES. 

C 

C HAKE  AN  OROERO  LIST  OF  PIVOTS. 

00  510  IslyN 
IN(I,2)=<1M(X.Z) 

IR=>IH(Iy4l 
510  lH<IR,3t=I 
G=G^GG 
GOTO  Z20 


THE  fOLLOHING.INSTKUCTlONS  IMPLEMENT  THE  FAILURE  EXITS. 

610  IF  (FIL.GT.O)  MRI TE f FI Ly 620 ) KylyJ 

620  FORMAT (//34XyZHELEMENT,17,10H  IS  IN  R0NtI5,llH  AND  C0LUNN,I5) 
G=-l 

GOTO  700 

630  IF  (FIL.GT.OI  MRITEI FILy 640)  I 

640  FORMAT (7/34Xy5HROM  «I5tl6H  HAS  NO  ELEMENTS) 

G=-2 

GOTO  700 

650  IF  (FIL.GT.O)  MRI TE ( F I Ly 660)  IRyJ 

660  FORMAT (//34Xt 35MTHERE  IS  MORE  THAN  ONE  ENTRY  IN  RON.I5y 
• IIH  ANO  COLUMN yI5) 

6=-3 
GOTO  700 

670  IF  (FIL.GT.O)  MRITE(FILy680) 

680  FORMAT (//34X,16H1A2  IS  TOO  SMALL) 

Gs-4 
GOTO  700 

682  IF  (FIL.GT.O)  MRITE(FILy 683) 

683  F0KNAT(//34X.16HIA1  IS  TOO  SMALL) 


..  BtSiriVAIUBlE  COPY 
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700  IF  (FIL.GT.O)  NMTeiFIL*  710) 

710  FOONAr(33N«ERftOR  kCTURN  FMON  INOANL  BECAUSE) 
720  RETURN 
ENB 
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SUBROUTINE  INOOPR  (A, INOZtIAZfN, IM,IP,0»B»GI 
IHPLICIT  REALMS  ( A-H , S>ZI • INTEGER  (I*RI 
IMPLICIT  REAL  t A-H, S -Z» «1NTEG£R  (I-RI 

INTEGER  IP  (N,  21, KP  (2),  10.(2) 

INTEGER»2  IN02 ( IA2I , IN  (N,5) 

INTEGER  IN02(1A*2),IN(N,5) 

DOUBLE  PRECISION  A(I A2 ), DIN) , 8 IN) 

DIMENSION  A<IA2),0(N) ,8(N} 

INOOPR  PERFORMS  THE  SOLUTION  OF  AX  > 3 NHEN  A NAS  BEEN  THROUGH  INOANL* 

IF  (G.LT.OI  GOTO  250 
DO  130  IlPsl,M 
IC1=1H(IIP,3> 

PP=IM(IC1,2I 

IF  (PP.EQ.O)  SOTO  130 

KPI1)=IP(IC1,1) 

KL(1>=KP(1)»IN(IC1,1)  -1 
B1=B(IC1> 

IF  (PP.EQ.l)  GOTO  100 
IC2=1NII1P*1,3) 

KP(2)=IP(1C2,1) 

KL(2)=KP(2)»IM(IC2,1)-1 

B2=B(IC2> 

IHIIC2,2>:0 
A12  = A(KP(  D) 

KP(1)=KP(1>»1 
B(IC1)=0(IC1) »81«A12*62 
8(IC2)=A12«31»0IIC2)*B2 
IQO  00  120  LL=1,PP 
KPP=KP(LL) 

KLL=KL(LL) 

IC-ICl 

IF  (LL.EQ.2)  lexica 
IF  (KPP.GT.KLL)  GOTO  120 
00  110  K=KPP.KLL 
IR=IN02IK) 

110  B(IR)=6<1R)-A(X)*8(ICI 

120  CONTINUE 
130  CONTINUE 

NRITEI6,1000)KLL 
1000  FO»MAT(1K,5H*****,IIO» 

00  200  IPIsl,N 

IIPsN*l-IPI 

IR2-IHIIIP,3> 

PP=IM(  IR2,2) 

BIRsO. 

KPPsIP(IR2,l) 

KLLsKPP*INIIR2,l)-l 
IF  (PP.EQV2)  KPPaKPP*! 

IF  (KPP.GT.KLL)  GOTO  150 
DO  149  K=KPP,KLL 
IC‘IN02(K) 

149  BlRsBIR-A(K)»a(IC) 

150  IF  (PP.ta.l)  3( IR2)«B(  IR2>/0(1R2) »dCR 
IF  (PP.EQ.O)  GOTO  160 

IF  (PP.NE.2)  GOTO  200 
A12sA(KPP-l) 
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lRJ«XM(XXP*i*Sf 
BCXIl2»«B(XltCI  »BaR2l*BlR 
BlXII3>>BIXRS)*B&2*aXR 
GOTO  200 

160  IRl>XN(XXP*tt3> 

B(1R2I»8(XR2)  *0tIR2>*BIR 
BCXRl)sB(XRll*0 (XPCXRl *11 I^BXR 
200  COMTXNUE 
250  RETURN 
ENO 
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SUBROUTINE  CONPRE ( At IRNt lA«Nt IM»IPtROMI 

this  subroutine  is  identical  to  the  HARRELL  SUBROUTINE  LAOSEO. 
DOUBLE  PRECISION  A(IA) 

OINENSION  A(IA) 

INTEGER»i»  IP(N) 

INTEGER  IP(N) 

INTEGER'2  IRNCI Al tlH(N) 

INTEGER  IRNCI A} tlN(N) 

LOGICAL  ROM 

CONNON/INOEFI/FIL.NUELtLRONtLCOL  tNCPtNUCL 
NCPsNCP*! 

00  5 

STORE  LAST  ELEMENT  OF  ENTRY  IN  lHiJt,  THEN  OVERNRXTE  IT  BY  -J 
NZ=IH( J) 

IF  iNZ.LE.B)  GOTO  S 
KsIP( J)»NZ>1 
IN(J)sIRN(K) 

IRN(K)s-J 
5 CONTINUE 

KN  IS  POSITION  OF  NEXT  ENTRY  IN  COMPRESSED  FILE. 

KNsQ 
IPI^O 
KL=NUCL 

IF  <ROH)  KL*NUEL 
C 

C LOOP  THROUGH  OLO  FILE  SKIPPING  ZERO  ELEMENTS  AND  MOVING  CENNINE 
C ELEMENTS  FORMARO.  THE  ENTRY  NUMBER  BECOMES  KNOMN  ONLY  MHEN 
C ITS  END  IS  DETECTED  BYPRESENCE  OF  A NEGATIVE  INTEGER. 

00  25  KsltKL 

IF  <IRN(K) .EQ.B)  GOTO  25 
KN*KN*1 

IFIRONI  A(KN)3A(K> 

IF  aRN(K).GE.a)  GOTO  29 
C 

C ENO  OF  ENTRY.  RESTORE  IRNIK),  SET  POINTERS  TO  START  OF  ENTRY  ANO 
C STORE  CURRENT  KN  IN  IPI  READY  FOR  USE  NMEN  NEXT  LAST  ENTRY  IS 
C DETECTED. 

J=-IRNIK) 

IRNfK|3lM(J> 

IP( J)=IPI»1 

IM(J)«KN-IPI 

IPI*KN 

20  IRN(KN)xIRN(K) 

25  CONTINUE 

IF  (ROM)  GOTO  26 

NUCLxKN 

GOTO  2T 

26  NUELsKN 

27  RETURN 
ENB 
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BLOCK  OATA 

CONNON/INOCFX/FXL »NUEL «LROM,LCOL »NCP»NUCL 
INrESEB  FIL 
OATA  FXL/B/ 

ENB 


REFERENCES 


1.  "The  NASTRAN  Theoretical  Manual,"  Edited  by 

R.H.  MacNeal,  National  Aeronautics  and  Space  Administration, 
Washington,  D.C.,  (1969),  p.  3.1-1. 

2.  Eisenstat,  S.C.  et  al.,  "Yale  Sparse  Matrix  Package 

I.  The  Symmetric  Codes,"  Yale  University  Dept,  of  Computer 
Science  Report  No.  112,  New  Haven,  Connecticut. 

3.  Eisenstat,  S.C.  et  al.,  "Yale  Sparse  Matrix  Package 

II.  The  Nonsymmetric  Codes,"  Yale  University  Dept,  of  Computer 
Science  Report  No.  114,  New  Haven,  Connecticut. 

4.  Munksgaard,  N.,  "Fortran  Subroutines  for  Direct 
Solution  of  Sets  of  Sparse  and  Symmetric  Linear  Equations," 
Institute  for  Numerical  Analysis  Rept.  No.  77-05,  The  Technical 
Univ.  of  Denmark,  2800  Lyngby,  Denmark  (Apr  1977) . 

5.  Duff,  I.S.  et  al.,  "Direct  Solution  of  Sets  of 
Linear  Equations  Whose  Matrix  Is  Sparse,  Symmetric,  and 
Indefinite,"  Atomic  Energy  Research  Establishment  Report  CSS  44, 
Harwell,  Dideot,  Oxon,  United  Kingdon  (Jan  1977) . 

6.  Rheinboldt,  W.C.  and  C.K.  Mesztenyi,  "Programs  for  the 
Solution  of  Large  Sparse  Matrix  Problems  Based  on  the  Arc- 
Graph  Structure,"  University  of  Maryland  Computer  Science 
Center  Technical  Rept.  TR-262,  College  Park,  MD  (Sep  1973). 

7.  Mesztenyi,  C.K.  and  W.C.  Rheinboldt,  "Further 
Programs  for  the  Solution  of  Large  Sparse  Systems  of  Linear 
Equations,"  University  of  Maryland  Computer  Science  Center 
Technical  Rept.  TR-394,  College  Park,  MD  (Aug  1975). 

8.  Gignac,  D.A. , "CSKYDG2:  An  Out-of-Core  Cholesky 
Algorithm  Equation  Solver  (with  Respect  to  Profile)  for  Large 
Positive  Definite  Systems  of  Linear  Equations,"  Naval  Ship 
Research  and  Development  Center  Rept.  4655  (Mar  1975) . 


35 


I 


